oo 

o 

O 
(N 

U 

Or 
< 



J2 



X 



Orbital-free effective embedding potential at 

nuclear cusp 

Juan Maria Garcia Lastra a , Jakub W. Kaminski, 

and Tomasz A. Wesolowski 

Universite de Geneve, 

Departement de Chimie Physique 



<D ; 30, quai Ernest- Ansermet, 

CH-1211 Geneve 4, Switzerland 
^ ! a Universidad del Pais Vasco 

Departamento de Fisica de Materiales 
j* • E-20018 Donostia-San Sebastian, Spain 

April 3, 2008 

> 

f— s | Abstract 

\Q ■ A new approach to approximate the kinetic-energy-functional de- 

pendent component (vt[pA, Pb\{t)) of the effective potential in one- 
electron equations for orbitals embedded in a frozen density environ- 
ment (Eqs. 20-21 in [Wesolowski and Warshel, J. Phys. Chem. 97, 

Q (1993) 8050]) is proposed. The exact limit for vt at pa — ► and 

f PBdr = 2 is enforced. The significance of this limit is analysed for- 
mally and numerically for model systems including a numerically solv- 
able model and real cases where f pBdr = 2. A simple approximation 
to vt[pA, Pb](t) is constructed which enforces the considered limit near 
nuclei in the environment. Numerical examples are provided to illus- 
trate the numerical significance of the considered limit for real systems 
- intermolecular complexes comprising, non-polar, polar, charged con- 
stituents. Imposing the limit improves significantly the quality of the 
approximation to Vt\pA, Pb](j*) f° r systems comprising charged com- 
ponents. For complexes comprising neutral molecules or atoms the 
improvement occurs as well but it is numerically insignificant. 



1 Introduction 

Numerical methods to study electronic structure in condensed matter use 
mainly techniques developed for periodic systems. In many cases, however, 
methods developed for finite systems are also used. They are especially 
adequate for ionic solids, liquids, molecular crystals, clusters of molecules, 
for instance, to study features of the electronic structure which are local in 
character. In such a case, the electronic structure is modelled only in some 
well-defined region in space of direct relevance. The effect of the atoms out- 
side of this selected region (referred to as environment in this work) is taken 
into account by some embedding potential. Different strategies are applied in 
practice to represent the embedding potential. They differ in the choice of 
descriptor of the environment. The roughest approximation is to neglect the 
environment entirely. Such simplification is commonly used to study chemical 
bonding and reactivity in condensed phase if the solvent in which the reaction 
takes place is known to play a secondary role. Representing the environment 
(discrete or continuous, polarisable or not) by the electric field it generates, 
makes it possible to take into account the effect of the environment [U, [2] . 
Such classical treatment of the effect of the environment on the electronic 
structure is commonly used both in chemistry and in materials science (for 
review see Ref. [3])- The embedding potential in such methods is obviously 
orbital-free. It is, however, not exact because the quantum statistics nature 
of electrons is completely neglected. Taking into account the fermion nature 
of electrons might proceed by following a similar strategy as the one applied 
by Phillips and Kleinman in the construction of pseudopotentials in order to 
eliminate explicit treatment of core electrons [4j. For recent developments 
along these lines, see Ref. [5]. 

Using the following elements of the Hohenberg-Kohn-Sham formulation of 
density functional theory: Hohenberg-Kohn theorems [6] , a reference system 
of non-interacting electrons [7j, and the corresponding density functional of 
the kinetic energy (T s [p]) [8] in particular, leads to the embedding potential 
which is exact in the limit of exact functionals and orbital-free i.e. does not 
involve other descriptors of the environment than its electron density [H] . The 
pure-state non-interacting f-representable electron density p™ m , such that 
added to some arbitrarily chosen density associated with the environment 
(Pb) minimises the Hohenberg-Kohn energy functional for the whole system, 
can be obtained from the following one-electron equations (Eqs. 20-21 in 



Ref. 0): 
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where pa = 2 J2i A |0«| an d an d v^ CED [pA, Pb', t\ denotes a local potential 
which depends only on electron densities pa and p#. The label KSCED 
(Kohn-Sham Equations with Constrained Electron Density) is used here to 
indicate that the local potential differs from that in Kohn-Sham equations [7] 
for either the total system {v ks [pa + Pb',^\)) or the isolated subsystem A 
{v KS [pA' ) r\)). Also the one-electron functions ({<f>f}) obtained from Eq. [T] 
are not optimal orbitals in neither Kohn-Sham systems. Atomic units are 
applied in all formulas which are given for spin unpolarised systems. 

The total effective potential in Eq. [T]is the sum of the conventional Kohn- 
Sham effective potential vfA [pa + Pb\ r\ for the whole system evaluated for 
the electron density p = Pa + Pb and another local potential (v t [pA, Pb](^)) : 



ksced 

J eff 



[pA, Pb] r\ = v eff [p A + Pb;t\+ v t [p A , Pb] (r) 
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where v t [pA, Pb\(j) involves functional derivatives of the functional T s [p\: 



v t {r) = v t [p A , pB](r) 
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Note that no restriction is made concerning the overlap between pa and pb 
in real space. The potential v t [pA, PB](r) can be alternatively expressed as: 



v t [pA,pB}(r) 
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where T' t ,[pa,Pb] denotes the following difference: 



7T>a, Pb] = T s [p A + P b] - T s [p A ] - T s [p B ] 



(5) 



For the sake of the subsequent discussions, it is convenient to split the 
total effective potential vf£ CED [pa,Pb]^\ into two components: the Kohn- 
Sham effective potential for the isolated subsystem A {vf^ [pA^'r]), which is 



pB independent, and the remaining part representing the environment: 

v5 s f CED [pa, pb; A = «5? [pa; A + vZ CED [pa, pb; A (6) 

where 

V^ CED [PA,Pb;A = vZttf+fg^d? (7) 
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where E xc [p] denotes the Kohn-Sham functional of the exchange-correlation 
energy j7]. 

Orbital- free effective embedding potential given in Eq. [TJ and its v t [pa, Pb] (r) 
component in particula, are used in various types of multi-level numeri- 
cal simulations (for a review, see Ref. pjj] or Refs. [UJ [121 CCS OH CGS OS 
IT?] [181 [HI [20] for representative recent reports). Such simulations deal 
with condensed matter systems, for which the electronic features of a se- 
lected subsystem (subsystem A) are subject to detailed investigation whereas 
Pb is subject to additional simplifications. Other formal frameworks use 
also v t [pA, PB]if) such as: Cortona's formulation of density functional the- 
ory [21], where ps is not an assumed quantity but a result of fully variational 
calculations [2TJ [221 1211 [231 [25] or linear-response time-dependent density- 
functional-theory description of electronic excitations localised in embedded 
systems [2"o] 12"?] . Finally, the orbital-free effective embedding potential given 
in Eq. [7J, and its v t [pA, Pb]{?) component in particular, are used in combi- 
nation with traditional wave-function based methods by Carter and collab- 
orators (see for instance Ref. |28j). For the formal analysis of applicability 
of such a combination, see Ref. [29] which shows that the exact embedding 
potential in such a case always comprises the v t [pA, Pb\(j) component. 

In practical applications, v t [pA, Pb}(^) is not used but some analytic ex- 
pressions approximating this quantity (v t [pA, Ps](^*)) for the obvious sake of 
practical advantages. This replacement results in errors in all derived quanti- 
ties which will be referred to as Vt-induced errors. In each case, Vt[pA, Pb](^*) 
is obtained by using analytic form of an approximated functional T s [p] into 



Eqs. [3J1SJ We will refer to such v t [pA, Pb\{t) as decomposable because the 
analytic form of all relevant quantities: T s [p], T™ ad [pA, Pb], and v t [pA, PB\(f)i 
is available. If the form of the used T s [p] comprises only low-level gradient- 
expansion [30] contributions, the corresponding decomposable Vt[pA, Ps](r) 
violates the exact limit for v t [pA, PB](r) at pa — > and / psdr = 2 (see 
Appendix A). Our interest in the local behaviour of v t [pA, Pb}^) at this limit 
is motivated by the fact that the corresponding conditionsoccur if p# com- 
prises two electrons tightly bound to a distant nucleus in the environment 
such as in the case of the helium atom, Li + cation, Be 2+ , etc. We expect 
that they are relevant also for heavier nuclei if in a volume element centred 
on the nucleus ps is dominated by a doubly-occupied orbital. 

The present work focuses on the investigation whether the considered 
exact limit is of any practical relevance. To this end, we apply the following 
strategy: i) We use a model system (Appendix B), for which the conditions 
Pa — ► and / psdr = 2 apply rigorously, to analyse the importance of 
enforcing the correct local behaviour of v t [pAi Pb\{^)- w) We construct a 
simple approximation to v t [pA, Pb\{t) obeying the considered exact limit in 
the vicinity of nuclei and analyse the numerical significance of imposing the 
considered condition in real systems where the conditions pa — > and 
J psdr = 2 do not apply rigorously. 

Our ultimate goal is a new approximation to v t [pA, Pb](^) which can be 
inexpensively evaluated in practice and obeys as much as possible of the 
relevant exact properties. It should pointed out in this context that the 
posit ion- dependency of Vt{pA, Pb\{t) is the result of non-homogeneity of pa 
and/or p#. Therefore, the symbol v t [pA, PB](r) (or Vt[pAi Pb]{t) if approxi- 
mated) is used throughout this work to indicate that this local quantity is a 
functional of pA and ps- Explicit position dependence is strongly undesired in 
density-functional-theory based methods because it is not straightforward to 
obtain i) such potential as a functional derivative of some density functional, 
and ii) functional derivatives of such potential needed in some formal frame- 
works [261 [27]. General symbols such v t or v t {f) are used in some discussions 
where the issue of explicit position-dependence is not relevant. 



2 Conventional (decomposable) strategies to 
approximate Vt\pA, Pb\{t) 

Before proceeding to the construction of the desired approximation to v t [p A , Pb] (t) 
obeying the considered limit, we overview the conventional construction of 
approximation to v t [pA, Pb](^*) and the local behaviour near a nucleus of the 
obtained potential. The conventional strategy, which is applied in our own 
works and the works by others so far, is to start from some explicit den- 
sity functional T s [p] and to use its analytic form to derive the corresponding 
approximate expression for T™ ad [p A , Pb]'- 

rr d [ PA , PB ] « fr d [ PA , PB ] = f s [ PA + PB ] - f s [ PA ] - f s [ PB ] (8) 

and to use the obtained analytic expression to obtain v t [p A , PB\(f) by means 
of functional differentiation. 
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v t [p A ,p B }(r) « v t [p A , p B }{r) = ~ s 
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(9) 
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This strategy can be applied for any approximated functional T s [p] provided 
its form makes it possible to obtain the analytic expression for v t [p A , /^(r). 
Simple functionals T s [p], which depend explicitly on densities and their gra- 
dients, are of particular practical interest. They lead to v t [p A , p B ](r) which 
depends explicitly only on p A and p B and their first- and second deriva- 
tives. In the original work by Cortona [21J, where the subsystem formulation 
of density functional theory was introduced, a decomposable v t [p A , PB\if) 
derived from the Thomas-Fermi [31] kinetic energy functional was used to 
study ionic solids. In our own works, only decomposable v t [p A) PB\(f) derived 
from gradient- dependent approximations to T s [p] were considered so far (see 
for instance the analyses of their accuracy in Ref. [32, 33] or their recent 
applications in multi-level computer simulations of condensed matter [T6]). 
Thomas- Fermi kinetic energy functional [31], which is exact for the uni- 
form electron gas, leads to the following approximate expression for T™ ad [p A , p B ] 

fr d ^[ P A,PB} = C TF J(( PA + PB f 3 - pf - pf) dr (10) 

where C TF = ^(3vr 2 ) 2 / 3 . 

The associated expression for vf F [p A , p B ]{r) reads: 

vf F [ P A, PB m = \c TF ((p A + PB f 3 - pf ) (11) 



Approximating T s [p] by means of the gradient expansion of the kinetic 
energy [3U] truncated to the second order leads to the following approximate 
expression for T™ ad [pA, Pb] [9]: 

f^GBAQ ^ pB] = T na d{TF) ^ p B ] _ J_ / ^JLZ ** Yff & (12) 

'^ J PaPb\Pa + Pb) 

The associated expression for v GEA2 [pa, Pb]{?) is given in Ref. [9]. 

For the the group of gradient-dependent approximations to T s [p] of the 
generalised gradient approximation form (32J [M] the analytic expression for 
T^IpatPb] reads: 

fr d{GGA) [pA,PB}=C TF J [(p A + Pb?" F(sab) - pTf(sa) - pTH*b)] df\13) 

where F GGj4 (s) (enhancement factor) depends on a dimensionless quantity 
s = 2(371-2^1/3 4/3 (reduced density gradient). Various analytic forms of F GGA (s) 
were proposed in the literature [351 ESI EU [38]. The associated analytic 
expression for vf \pAi Pb[(t) is given in Ref. [23] . It is worthwhile to notice 
that the GGA form is flexible and includes T TF [p] and T GEA2 [p] as special 
cases. Numerical values of s provide useful information about shell structure 
and the distance from the nucleus in atoms [39]. For an atom, s is known to 
be small near the nucleus, reach the values of about 3 in the valence region, 
and diverge exponentially to +oo at large distances. In molecules, it behaves 
similarly with a noticeable exception of stationary points of electron density 
(bond midpoints for instance) where s = 0. Each approximated functional 
given in Eqs. [T0iri3l comprises a dominant Thomas-Fermi component and 
satisfies two exact conditions: 

• T s [p A + Pb] - T s [p A ] - T s [p B ] = for non-overlapping p A and p B . 

• For uniform pa and ps, they recover the exact analytical expression for 

T s [p A + Pb] - T s [p A ] - T s [p B }. 

The common feature of each among the above approximations for T s [p] is 
that none of them yields the exact analytic form of v t [pA, Pb\{x) & t Pa — ► 
and / psdf= 2 (see Appendix A): 

V t [pA,PBW)-^V l r t [ P BW) = il^C - \^ (14) 

X Pb 4 Pb 



The expression given in Eq. [TTJ which provides the dominant contribu- 
tion to gradient-expansion based approximations to v t [pA, PB\if) does not 
comprise the relevant term at all whereas the second-order term provides 
only 1/9 of the exact expression. Figure [I] shows v emb [Pa, PB](r) for 

a spherically symmetric case: ps = PHe, v f x t(^) = — 2/r, and pa — ► 0, 
which represents a helium atom far from subsystem A. In the Figure as well 
as in the following discussion, r denotes the distance from the considered 
nucleus. The potential in the figure shows the features which are common 
also for heavier atoms if gradient expansion based approximations to T s [p] 
are used to derive v t [pA, PB\{r)'- a very narrow and deep well (reaching — oo) 
centred on the nucleus and surrounding it repulsive shell. The v t \pAiPB\if) 
component of the shown embedding potential is finite at the nucleus instead 
of behaving as - . Note that the exact term has the same form as the poten- 
tial due to Coulomb attraction by the nucleus of the charge Z. Therefore, 
the exact v t partially compensates this attraction to some extend because 
( is smaller than Z [40] . For the particular case considered in the Figure, 
the missing - component does not lead to any bound states (Appendix B). 
In general, however, the wrong asymptotic of the singularity at the nucleus, 
can lead to an unphysical transfer of electron density from the investigated 
system to its environment (charge-leak [2]). This can occur if the artificially 
attractive (not sufficiently repulsive) approximation to the orbital-free effec- 
tive embedding potential generates a bound state in the environment of the 
energy which is lower than the eigenvalue of the highest occupied embedded 
orbital associated with embedded subsystem A. Numerical cases confirming 
such scenario are known [121 03] . Moreover, the numerical solution of the 
Schrodinger equation with v emb [pa-, Ab](^*) for ps = Pu+ and the ex- 

ternal potential of the — - form, shows a deeply lying node-less bound state 
of the energy -0.209665 hartree with the maximum of the radial electron 
density at r max = 2.912 bohr (Appendix B). The above observations indicate 
that decomposable v t [pA, Ps](r) obtained from low-order gradient-expansion 
based approximations to T s [p\ might not be adequate for, at least, Li + cations 
in the environment. The fact that the bound state associated with an atom 
in the environment is too tightly bound to the nucleus indicate, however, 
that the problem might be also present in atoms comprising more electrons. 

In a subsequent section, a simple approximation to v t [pA, PB\(f) is con- 
structed based on these observations. In principle, the exact limit should 
be applied in any volume element in which pa vanishes and ps is obtained 



from a doubly occupied orbital. In the proposed construction, the exact 
limit is imposed only at volume elements near heavier-than-hydrogen nuclei 
expecting that the considered condition is most relevant there. 

3 Building-in the exact limit for v t [pA, Ab](^) 
at pA — ► and J psdf= 2. 

The aforementioned flaws of decomposable strategy to construct gradient- 
and Laplacian dependent approximations to v t [pA, Pb](^) suggest a bottom- 
up approach in which v t [pA, Pb\{t) is directly a target. A given approximated 
potential v t [pA, Pb](j) will be referred to as non- decomposable if the analytic 

form of its two individual components . and . cannot be 

Sp p= PA +p B s P p=PA 

reconstructed. The non-decomposable strategy is motivated by the fact that 
there are exact properties of Vt[pA, Pb](^*), which can be taken into account 
quite easily in v t [pA, Ps](^), whereas building-in them into some approximate 
functional T s [p] is less straightforward. Abandoning the decomposable strat- 
egy is motivated also by the results of our recent dedicated studies of the ac- 
curacy of various gradient-dependent approximations to T™ ad [pA, Pb], which 
revealed that there is no correlation between the accuracy of T™ ad [pA, Pb], 
v t [pA, PB\(f) and the errors in the parent gradient-dependent approximation 
to T s [p] [33]. It should be also pointed out that, that the individual contri- 
butions T™ ad [pA, Pb] are not needed in practice. 

The non-decomposable approximation to v t [pA, Pb](t) is constructed by 
enforcing the following exact conditions into its analytic form: 

• f™ d [p A ,p B ] — ► fr d{LDA) [pA,pB] for uniform p A and p B . 

• T™ ad [pA, Pb] — ► for non-overlapping p A and p B . 

• Vt[pA,pB] — ► v l t imit [p B } &t p A — ► and J p B df= 2. 

The first two conditionsare automatically satisfied by the decomposable 
gradient-expansion based approximations discussed in the previous section. 
Since such approximations proved to be sufficiently accurate for many sys- 
tems the same conditions are retained in the new construction. The last 
condition is the key element of the present construction. 

Before proceeding to the construction of the approximation obeying the 
considered exact condition we note that v t [pA, Pb](^) can be alternatively 



expressed as: 

VtiPA, Pb] = vt ecomposable [p A , p B ] + f[p A , p B ] ■ v l r it [ PB ] (15) 

All functionals in the above equation are determined locally and the argu- 
ment f is not written explicitly for simplicity and the functionals Vt[pA,PB\ 

r 1 -decomposable r -i 

and f[p A ,p B } are simply related (f[p A ,p B ] = ^'^"fr^,^] [pa ' Pb] if 

v l t imit [pA, Pb] is non-zero). The above form of v t [pA, Pb] provides a convenient 
for construction of approximation. It can be used for any decomposable ap- 
proximation to v t [pA,PB], which violaties the considered condition, and the 
functional f[pA,PB] has a clear physical meaning as a switching factor de- 
termining whether it is needed to add locally the missing component of the 
embedding potential. As far as the choice for the decomposable component, 
both the gradient-free potential given in Eq. [TTJand the decomposable poten- 
tial derived from the Lembarki-Chermette [36] approximation to T s [p] were 
shown in dedicated studies [221 S3] to be reasonably accurate if pA and p B do 
not overlap strongly. Both these approximations comprise the zeroth order 
contribution. If vf F [pA, Pb] is used as decomposable component in Eq. [151 
the terms vl tmtt [p B ] and y t ecom P° sa e [p A) p B ] are local functions depending ex- 
plicitely on p A , p B , Vp#, and V 2 p B . Approximating /[pa^Pb] by a local 
function depending explicitely on these quantities applied in Eq. [15] leads to 
an approximated potential requiring a similar computational effort as con- 
ventional low-order gradient-expansion based decomposable approximations 
to v t [pAiPB\- The first approximation made here is replacing the switching 
factor defined in Eq. [15] by a switching function: 

f[pA,pB}(r) ~ f(p A ,pB,Vp B} V 2 p B ) (16) 

The above considerations lead to the following general form of the ap- 
proximation to v t [pA, Pb]{?)'- 

v t [PA, Pb] = vJ F [p A , Pb] + f(p A , Pb, Vp B , V 2 p B ) ■ v l t imit [p B ] (17) 

The above general form provides a clear interpretation for the switching 
function, which can be used as guideline in construction of approximations 
- it "detects" such volume elements for which the conditions pa — > and 
/ p B dr = 2 are most relevant. 
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3.1 The switching function / for environments com- 
prising one-nucleus and two-electrons 

In constructing f(pA,PB, Vp^, V 2 ps) the following additional requirements 
(simplifications) are made: 

• f(pA, Pb, Vps, V 2 ps) is one in the vicinity of a nucleus to account fully 
for the missing - component. 

• The criterion for determining the range at which v l ^ mit [pB] is nuclear 
number independent. 

• f(pA, Pb, Vp#, V 2 /9b) does not depend on p& (to obtain the analytic 
form of fr d [p A , Pb]: (f[p A , Pb] « f\p B \) 

The above criteria are very restrictive and leave us with not many choices. 
The last one leads to the following form of the switching factor: 

f[pA,PB]-f(pB^p B ,V 2 PB ) (18) 

Approximating J[pa,Pb] by some function f(ps) is one of possible further 
simplifications. It is, however, very unlikely that a p^-based switching func- 
tion could be universal. The electron density near the nucleus depends on 
the effective nuclear charge ( [ID] and varies strongly from atom to atom. It 
is possible, however, to design an (^-independent criterion. To this end, we 
consider the reduced density gradient (s#(f)) defined as: 

IVpbI , 1fV . 

SB = 77T (19) 

2(37r 2 )V3 p ;f 

For pb obtained from hydrogenic orbital Is defined by some effective nuclear 
charge (, the following (^-independent observations can be made: i) For p 1 ^ = 
2|ls| 2 , s r=0 = (67r)~ ' = 0.376, and ii) v l ^ m%t [p 1 ^] changes sign from positive 
to negative at s r=2 /( = exp(4/3) • s r=0 = 1.426. These observations suggest 
that the switching function can take a very simple and ^-independent form 

Kpb, Vp B ) = f(s B ) = e(s B - sT n ) x e(sT x - a B ) (20) 

where Q(x) = 1 for x > and Q(x) = for x < and s™ n = 0.376 and 

g rnax = 1A2 Q. 
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Since Eq. [T] can be used also to obtain forces [25] , it is preferable to use 
a smooth switching from to 1 instead of G in the above definition. The 
simplest form of such a switching functional has the Fermi-Dirac statistics 
form: 



-l 



/ = (exp(A(-s B + st n )) + l) x (l - (exp(A(- SB + s^)) + 1)~ 

(21) 

where the parameter A determines the smoothness of the switch. 

A = 500 in Eq. [21] leads to equivalent results to that obtained with the 
step function O(x) (differences in dipole in the range of 10~ 6 Debye and 
orbital energies in the range of 10 -10 hartree). All results discussed in this 
work are obtained with A = 500. Smaller values corresponding to even 
"softer" switching can be also used in our numerical implementation but is 
not considered here in order to minimise the use of adjustable parameters. 

A switching factor / constructed following the above restrictions can be 
used to investigate the importance of enforcing the exact limit for v t [pa, Pb\ if) 
at pA — > and / psdr = 2 but only in cases where the environment com- 
prises one nucleus and two electrons. Not only the condition / psdr = 2 
apply rigorously but the considerations leading to the (^-independent values 
of s t q ix and s B wn apply as well. Such model systems as Li + -H20 and Be 2+ - 
H 2 complexes if the water molecule is considered as subsystem A and the 
cation as subsystem B (environment) fall into this category. At equilibrium 
geometry addition of / ■ v t tm lt [pB] results in a desired effect on calculated 
properties for these complexes. The lowest unoccupied embedded orbital as- 
sociated with subsystem A is indeed localised on the cation and its energy is 
shifted by 0.262 eV in the case of Li + and by 0.830 eV in the case of Be 2+ . 
Addition of / ■ v l ^ m%t \pB\ reduces also the dipole moment of a water molecule 
in the vicinity of the cation by 0.071 Debye and 0.409 Debye for Li + and 
Be 2+ , respectively. Such noticeable numerical effects obtained for systems, 
for which the condition / psdr = 2 applies rigorously, indicate clearly that 
the exact limit considered might be relevant for practical calculations where 
the environment is larger. 
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3.2 Fine-tuning of the thresholds in the switching func- 
tion / 

Using only universal parameters s™ n and s^^ is very appealing but the 
reasoning leading to their numerical values does not apply in real systems. 
It allows one to study the importance of the considered exact conditions 
but only in particular cases. The approximation to v t [pA, Pb}(^) defined in 
Eq. [T7] and using the switching function given in Eq. [21] has, therefore, little 
practical value. The construction of the switching function described in the 
previous section and choice of the thresholds s™" and s^ ax in particular can 
be expected not to be adequate for other systems for the following reasons: 

• The criterion 0.376 < sb < 1.426 applies at hydrogen nucleus where 
addition of ff™'[ps] is not expected to be needed. We recall here that 
it is the danger of a collapse of electron density on a doubly occupied 
hydrogenic Is orbital provides the physical motivation for introducing 
the v l t imit [pB]- Moreover, numerical studies on molecular electron densi- 
ties indicate clearly that enforcing the local behaviour of the density of 
the kinetic energy near nucleus corresponding to the von Weizsacker ex- 
pression leads indeed to significant improvements of the approximation 
to T s [p] for all nuclei except that of hydrogen [4T] . 

• For heavier atoms, electron density at the nucleus comprises contribu- 
tions from other orbitals than the hydrogenic Is. This can lead to the 
possibility that is v l t zmit [pB] is negative although sb < 1.426. 

• The criterion 0.376 < sb < 1.426 might also be satisfied near stationary 
points of the electron density such as bond midpoints. It is very unlikely 
that the condition / psdf = 2 can be relevant to any volume element 
centred on a stationary points. Therefore addition of v l ^ m%t \pB\ lacks 
formal justification there. Although the condition 0.376 < s# assures 
that v^ mit \pB\ is not added at the stationary point, where |V/?b| = 0, 
or in the close proximity to it, the numerical value of this threshold 
requires verification in real systems.lt should be added at this point 
also that, if linear combination of atomic orbitals is used to construct 
embedded orbitals, the quality of description of the density at the nu- 
clear cusp depends on the used basis set. A weaker criterion should 
be used in practice to assure that the v l t zmit [pB] is indeed added in the 
vicinity of the nucleus. 

13 



We start with the choice made for s^ ax . Numerical analyses in the sys- 
tems discussed in the next section show that v l t imit [pB] is negative locally even 
if sb is smaller than 1.426 for heavier nuclei. This suggest that this thresh- 
old should be reduced. The smaller is the value for this threshold the less 
probably is inclusion of negative v^ TOii [ps], which is desired from the point of 
view of universality of this threshold, but it comes at the expense of loosing 
a part of the desired effect at two-electron nuclei by reducing the range at 
which addition of ff mit [ps] applies. We use the model system considered 
in Appendix B to estimate the effect associated with the reduction of this 
range. The underlying assumption leading to the value of 1.426 value is rig- 
orously true in the model system. The desired effect of reducing the charge 
distribution on top of the nucleus is achieved mainly by adding v l t irntt [pB](r) 
very close to the nucleus i.e. where sb < 0.6. Increasing further the value of 
the sb threshold leads to smaller effect. Moreover, adding locally v l t imit [pB] 
to the potential near the nucleus leads to negligible effect on orbital energies 
in the system considered in Appendix B (it reaches a peak of about 10~ 4 eV 
at sb = 1.426). These results for the model system indicate that any choice 
for 0.6 < s^ ax < 1.426 is acceptable. In real systems discussed in the next 
section, lowering the threshold from 1.426 to to 0.9 assures that v l t imit [pB] is 
added only if it is positive. 

The fine-tuning of s B vin follows other considerations. We note that, in the 
model system considered in Appendix B, s™ n can be reduced even to zero 
without affecting the results because lower values of sb than 0.376 do not 
occur near the nucleus. To make sure that no nuclei is overlooked even if 
the chosen atomic basis set in practical calculations is such that the exact 
relation sb = 0.376 for r — ► cannot be rigorously satisfied, s™ n is reduced 
from 0.376 to 0.3. This change leads to negligible numerical effects if ps 
corresponds to atomic electron densities. For molecular ps, retaining the 
criterion based on the s™ n is necessary to avoid unjustified additions of 
v l t imit [pB] near stationary points. 

The criteria based only on s B im and s^ ax are not sufficient if the environ- 
ment comprises hydrogen atoms because they are satisfied also at hydrogen 
nucleus. To avoid adding the v l t imit [pB] near hydrogens, the proposed switch- 
ing function includes additionally the criterion based on smallness of ps- It 
is required that pb is larger than the square of the Is wave function of the 
hydrogen atom (Z=l) at r = which equals to 1/ir = 0.318. Concerning 
p™ n , increasing the idealised value of 0.318 to even 1 does not affect the 
results for hydrogen-free systems because the density on top of any nucleus, 
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which is heavier than hydrogen, is at least one order of magnitude larger. 
Increasing the value p™ n is desired for the same reasons as the ones moti- 
vating the decrease of s B nn . The value of p B vin = 0.7 was arbitrary chosen for 
practical calculations. 

The final form of the "fine-tuned" switching function of more general 
applicability and used in the subsequent section for studying the importance 
of imposing the considered exact limit in real systems take the following form: 



/ = (exp(X(-s B + s^ m )) + iy x(l-(exp(X(-s B + s r r)) + iy 

x (expM-PB + pf* )) + l) _1 (22) 

where s^ in = 0.3, s^ ax = 0.9, p^ in = 0.7. 

Eqs. [TTj [TH [TTJ, and [22] define the potential which will be referred to as 
v^ dsd [pa, Pb], (Non-Decomposable approximation using first- and Second 
Derivatives of p). The notion of non-decomposability is brought up here be- 
cause the the second term in Eq. [T7] does not have the form of a difference 
between two functional derivatives of some common explicit density func- 
tional T s [p\. The analytic expression for T™ ad ( NDSD ^{pA, Pb], which yields 
v^ dsd [pa, Pb] after functional differentiation with respect to pa can be eas- 
ily constructed. Its decomposable component is given in Eq. [10] and the 
non-decomposable v l t zmit [p B ] component is pa independent. The functional 
generating vf DSD [pA, Pb] reads therefore: 

fr d{NDSD) [PA, PB ] = c TF J((pA+p B f 3 -pT-pT)df (23) 

+ J f(p B , V PB ) ■ PA(?)v l r t ip B }(r)dr + C[p B ] 

where C[p B ] is p^-independent. To assure the proper dissociation limit C[p B ] 
must vanish. 

4 Numerical validations 

Procedure to analyse {^-generated errors 

In practical applications of Eq. [T] the results depend on p B as well as on the 
used approximation to v t [pa, Pb] if) ■ As far as the quality of the used approxi- 
mation to v t [pA, PB]if) is concerned, a general procedure was proposed in one 
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of our earlier works [33J EI] • Its principal element is the comparison between 
numerical values of the calculated property (energy components, dipole mo- 
ments, total electron density, etc.) obtained from two fully variational formal 
frameworks: that of Cortona [21] and that of Kohn and Sham [TJ. Results 
obtained from both frameworks are not exact but the difference between 
them can be attributed only to the approximation used for v t [pA, Pb\{t) if ah 
technical parameters (approximation to the exchange-correlation functional, 
basis sets for expanding orbitals, algorithms to calculate used matrix ele- 
ments) are the same. We point out here, however, that direct comparisons 
between the total electron densities derived from Kohn-Sham- and Cortona's 
calculationsare cumbersome because these quantities are local. In practice, 
it is more convenient to use global quantities (norm of the difference be- 
tween these densities, or selected observables) in such analyses (see the next 
section). 

To obtain the pair of electron densities pa and pb which minimises the 
total energy in Cortona's type of calculations, a self-consistent super cycle of 
embedding calculations (freeze- and-thaw cycle) is performed. At each itera- 
tion, Eq. [TJ are solved. In the subsequent iteration, pa and ps exchange their 
role in Eqs. [TjThe freeze- and-thaw iterations continue until self-consistency 
In the end, a pair of electron densities (p A and p B ) and the corresponding 
two sets of embedded orbitals is obtained. Obviously, the notion of embedded 
system and its environment becomes meaningless because both subsystems 
are treated on the equal footing. Freeze- and- thaw calculations are conducted 
in practice for small model systems to validate the used ps in large scale 
multi-level numerical simulations [16] or, as it is made in the present work, 
to asses the used approximation to v t [pA, Ps](^*)- 

The {^-generated errors in the complexation induced 
dipole moments due to violation of the limit for v t at 

Pa — ► and / psdf '= 2. 

■^-generated errors in complexation induced dipole moments can be expected 
to be strongly affected by the local behaviour of the used Vt[PA, Pb\{t) near 
nuclei. Lack of sufficient repulsion near the nucleus might lead to an artificial 
transfer of electron density between subsystems reflected in the numerical 
values of the dipole moment. Therefore, this quantity was chosen for the 
analysis of the ^-generated errors in a representative set of intermolecular 



16 



complexes including charged, polar and non-polar ones at their equilibrium 
geometries. For key details of the numerical implementation of the relevant 
equations, see Ref. [4"9] . 

Tables [1] and [2] collect the complexation induced dipole moments in neu- 
tral or charged complexes, respectively. First of all, switching on the v l g mit \pB\ (r) 
term decreases the {^-generated errors in each of the considered cases. For 
systems comprising neutral subsystems the effect on the errors are negligible. 
This indicated that the origin of the errors lies not in the violation the con- 
dition considered in this work. For systems comprising charged components, 
the effect of imposing the considered limit is evident. The errors are invari- 
ably reduced. The reduction of the relative errors depends on the system 
from such a case as Li + -H 2 (from 9.8% to 7.9%) to a reduction by factor 
2 or 3 in Na+-Br~ and (from 1.4% to 0.7%) and Na + -H 2 (from 0.23% to 
0.07%). Similarly as in the group of complexes formed by neutral molecules, 
the origin for the remaining errors lies somewhere else. The above numerical 
examples lead to the following principal conclusions: 

• Violation of the limit for v t [pA, PB](r) at pa — ► and / psdr = 2 
contributes to the overall error in the calculated quantities but this 
contribution varies from one system to another. It is rather negligible 
for complexes formed by neutral components. It is numerically signifi- 
cant for complexes comprising charged components. 

• Our simple strategy to impose the considered limit locally in the vicin- 
ity of nuclear cusps leads invariably to reduction of errors. Therefore 
it can be used generally as correction to any approximation violating 
the above limit. 

• The construction of the approximation obeying the considered limit, 
and the used switching criteria in particular, corresponds to a real 
case where a distant nucleus is surrounded by a frozen- density shell 
comprising two electrons (He, Li + , Be 2+ , etc.). The fact that the errors 
are reduced also for systems, where these idealised conditions do not 
apply, indicates that the considered condition is important and should 
be taken into account in construction of approximations to v t [pA, Pb\{t)- 

• For systems, where £> t -generated errors were not reduced by imposing 
the considered exact limit, their origin must be looked for somewhere 
else. 
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The effect of imposing the limit for v t at pA — ► and 
/ psdr — 2 on orbital energies 

In this section, we analyse the complexation induced shifts of orbital energies 
derived using the approximated potential v^ dsd [pa, Pb\{t) considered in the 
previous section. Opposite to the dipole moments discussed previously, direct 
comparisons between the calculated shifts and the corresponding reference 
data are less straightforward. However, we investigate the numerical effect 
associated with imposing the exact limit for v t [pA, Pb\{j) at Pa — ► and 
/ psdr = 2 in view of the numerical practice which indicates that shifting 
the levels of unoccupied orbitals localised in the environment would be de- 
sirable. Conventional decomposable approximations to Vt\ftA, Pb](t) lead to 
artificially low levels of unoccupied orbitals in the environment which might 
cause unphysical effects such as charge-transfer between subsystems (13] or 
erroneous other observables |42j . 

The numerical results for a model system considered in Appendix B show 
that inclusion of the v t imit [pB] term everywhere where 0.6 < sb < 0.9 into 
the effective embedding potential leads to a positive shift of the energy level 
of the unoccupied orbital. In real intermolecular systems, the conditions 
considered in Appendix B (pa — >■ and ps = 2|ls| 2 ) are not satisfied. The 
subsystems are in finite separation and the use of atom-centred basis sets for 
each subsystem, which include all atoms (supermolecular expansion labelled 
as KSCED(s) in Ref. [IS]), results in the fact that pa can be significant at 
a nucleus associated with subsystem B. For the same reasons and the fact 
that the considered nuclei include also atoms with occupied 2s shell, also the 
second assumption is not satisfied rigorously. Therefore, it is useful to verify 
in practice to which extend inclusion of the v l t irntt [pB} term affects the orbital 
levels if these asymptotic conditions do not apply. 

Tables |3] and 0] collect the values of energy levels corresponding to the 
lowest lying orbital localised mainly in the environment for the previously 
considered complexes. The calculations are made not in the end of the freeze- 
and-thaw cycle but for the same ps obtained from Kohn-Sham calculations 
for the isolated subsystem B. Including v l t tmtt [pB] leads to the shifts of the 
energy levels of unoccupied orbitals of the magnitude which significantly 
larger than that in the model system. It is a very desired effect of the new 
approximation. 

In the Li + -H20 case discussed in Ref. [42J, the energy of the lowest unoc- 
cupied embedded orbital localised on Li + crosses that of the highest occupied 



embedded orbital localised on H 2 at the intermolecular distance of 13 A if 
Eq. [TT]is used for v t [pA, Pb](j)- At larger separations, the self-consistent pro- 
cedure to solve Eqs. CD does not converge due to localisation of the highest 
occupied embedded orbital which jumps between subsystems in subsequent 
iterations. Addition of the v l t imit [pB] term shifts the energy of the unoccupied 
embedded orbital localised at Li + . As the result, no crossing of levels occurs 
even at intermolecular separations as large as 18 A. The occupied levels are 
affected less strongly (see Tables [5] and [6]) . 

5 Discussion 

In principle, any decomposable approximation can be used as the first term 
of Eq. [T7] instead of the term derived from Thomas- Fermi functional. In 
the approximation introduced in this work, the decomposable component 
of vf DSD [pA, Pb] is gradient-free and does not contribute to the asymptotic 
local behaviour of v t [pA, Pb]^ at the nuclei. Only the second term en- 
forces the desired behaviour. Should a gradient- dependent alternative for 
the first term be considered, a proper care should be taken to avoid double- 
counting of v l t tmit [pB] at a nucleus. For instance, vf EA2 [pA, Pb] comprises 
already 1/9 of v l t imit [pB]- To verify whether further improvements are pos- 
sible following this lines, two approximations were considered by replacing 
vJ f [Pa,Pb] in Eq. El by either v? EA2 [p A , Pb] [9] or ^ GA97 [p A ,p B ] [32]. In- 
cluding the v l t tmtt [pB] term into v s [pa, Pb]{t) brings improvements in both 
cases (see Table [7]). This indicates v l t imtt [pB] should be enforced universally 
on any approximation to v s [pa,Pb]{?)- Compared to the results obtained 
using v^ dsd [pa, Pb]{^) discussed earlier in this work, the two alternative 
non-decomposable approximations v s [pa, Pb]('t) are not better. 

We recall here the main reasons for singling out v GGA97 [pa, Pb]) among 
other decomposable ones which depend explicitly on densities pa and ps as 
well as their first- and second derivatives. 

• vf EA2 [pAi Pb] obtained from the second-order gradient-expansion ap- 
proximation leads typically to worse results than that obtained from 
zeroth order [33], H5] indicating that the contribution to v t [pA, Pb] due 
to the second term in Eq. [T2] is erroneous. The deterioration of the 
results is pronounced the most at very small overlaps between pA and 
Pb- Note that in the present work this flaw of v GEA2 [pa, Pb] manifests 
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itself in the absence of convergent solutions of Eq. [T] in Na + -Cl and 
Li + -H 2 cases (see Table [7]). 

• v GGA97 [pa, Pb] was introduced as a pragmatic solution replacing vf EA2 [pa, Pb] 
Due to its analytic form, the gradient-dependent contribution disap- 
pears at small overlaps between pa and ps- 

• The functional T LC94 [p] generating the decomposable approximated po- 
tential vf GA97 [pA, Pb] is known to be a very good approximation to 
T s [p\. 

The above reasons for singling out v GGA97 [pa, Pb] are not applicable for the 
non-decomposable construction presented in this work. In v^ dsd [pa, Pb] the 
problematic second-order term lying at the origin of flaws of v GEA2 [pa, Pb] is 
either present where it is needed to assure the correct asymptotic limit (i.e. 
in the vicinity of nuclei) or absent. Numerical results collected in Tables HJ 
121 and [7] support fully the above formal reasons to consider vf DSD [pA, Pb] as 
the successor of v G GA97 [pa, Pb]- 

6 Conclusions 

Enforcing the considered exact limit for vt[pA, Pb]^), as it is made in the 
proposed approximation, leads to a significant reduction of {^-generated er- 
rors for charged systems. Errors in the complexation indiced dipole moments 
are reduced by more than 50% in some cases. For neutral systems, reduction 
of error takes also place but its magnitude is typically negligible. This indi- 
cates that one of important sources of inaccuracies in the conventional (i.e. 
decomposable and gradient-expansion based) approximations to v t [pA, Pb]{?) 
was identified. The origin of the remaining contributions Wj-generated errors 
lies probably somewhere else. The analytic form of the component of the lo- 
cal embedding potential enforcing the correct considered limit is, indeed, an 
approximation to v t \pA, Pb\{t) because its position dependency is indirect - 
through the density and its gradient. The function which was used to switch 
on the exact limit was designed based on the analysis of model system of 
relevance for elements of the first-, second- and the third period. Using the 
local potential introduced in this work (vf DSD [pA, Pb](?) given in Eq. [TTJ as 
an alternative to potentials derived using conventional strategy of deriving 
it from analytic form of functionals based on low-order terms in the gradient 
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expansion of T s [p] [30] is recommended for the following reasons: 

i) Formal: We believe that a proper strategy to improve approximations 

to functionals in density functional theory should proceed by imposing the 

most relevant exact conditions and v^ dsd [pa, Pb\{t) was constructed in this 

way. 

ii) Practical: The numerical results reported in this work show indeed that 

imposing this condition improves the obtained electron density and that the 

improvement varies from negligible to significant depending on the system. 

Hi) Numerical: Evaluating v^ D [pa, Pb\{x) involves the same quantities 

as evaluating its counterparts derived from gradient-expansion- (up to second 

order) and so called generalised gradient approximations to T s [p\. 

Concerning the area of applicability of v^ DSD \pAi Pb]{x)i it should be 
underlined that some arbitrary choices were made concerning the criteria for 
"detecting" the vicinity of a nucleus based only on electron density p# and its 
derivatives in the construction of this approximation. The chosen criteria are 
most adequate for such nuclei, where the total electron density is dominated 
by the hydrogenic Is orbital. 

The local potential vf DSD [pA, Pb]{t) introduced in this work is intrin- 
sically non-decomposable. Although the analytic form of the functional 
rpnad{NDSD) ^^ ^j anc j ^ e potential v^ dsd [pa, Pb] (r) are given in this work, 

the analytic form of neither T^ DSD [p] nor — s —g — — is available. Therefore, 
the introduced here non-decomposable strategy can be seen as the first at- 
tempt to decouple the search for T s [p] and its functional derivative, which 
are needed in orbital-free calculations, from the search for an adequate ap- 
proximation for the kinetic-kinetic-energy dependent component of the effec- 
tive orbital-free embedding potential. Opposite to orbital-free strategy [3T] . 
neither T s [p] nor its functional derivative are needed in methods applying 
orbital-free effective embedding potential given in Eq. [7J 

Finally, imposing the exact limit for v t [pA, Pb\(j) at pa — > and / psdf = 
2 leads to shifts of the level of the unoccupied orbitals localised in the en- 
vironment. Such shift is strongly desired in view of the earlier reports on 
possible practical inconveniences resulted from artificially low position of 
such levels when approximations not taking into account the cusp condition 
are used [321133]. 
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Appendix A 

Orbital-free embedding potential for pa — ► and p B = 2\ls\ 



For small Sp such that 5p — ► 0, 



5Tr d [ P A,p B ] = Tr d [p A +5p, PB ]-Tr d [PA,PB] 



(24) 



T s [p A + 5p + p B ] - T s [p A + p B ] - T s [p A + Sp] + T s [p A ] 



ST s [p] 
Sp 



Sp(r)df 



P=PA+PB 



ST.\p] 
Sp 



Sp(r)df+ 0{S 2 p) 



P=PA 



If also pa is small i.e. p A — ► 0, 



TT d [p A + s P , Pb ] - rr d [PA, pb] = J 



Sp 



Sp(r)df+ 0{S 2 p) 



P=Pb 



Therefore, 



5Tr d [p,p B } 

Sp 



p^O 



ST a \p] 
Sp 



(25) 



(26) 



P — >Pb 



The above result that the kinetic-energy component of vf^ ED is just the 
functional derivative of the T s [p] calculated for p = p B makes it possible to 
express it analytically for any ps which comprises just two electrons. For one- 
electron and two-electron- spin-compensated systems the exact expression 
reads |Sj: 



T 8 [p\ = TY\p] = J l -^^dr for J pdr = 2 



Therefore, 



g7M 

Sp 



P=PB 



l \Vp B \ 2 lV 2 p B 
8 P 2 B ~4 p B 



if / psdr = 2 



(27) 



(28) 
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Using Eq. [2S] in Eq. [2S] leads to the asymptotic form of the kinetic energy 
component of v^^ ED in the case where pa and p# do not overlap signifi- 
cantly and pb is a two-electron spin-less electron density reads: 



5Tr d iPiP B ] 



Sp 



_1 \V P b\ 2 lV 2 p B 

P=PA -^0j PB dr=2 8 Pb 4 Pb 



(29) 



Pb representing a doubly occupied hydrogenic Is function (Is = a/C 3 / 71 " ' 
exp(— (r)) reads: 

Pb (»0 = Pb (0 = 2 • C 8 /* • exp (-20") , (30) 

For p B s , Eq. [29] leads to the following potential: 
ST^ip, p B }_ {f) 

P=PA — >0,pb=Pb 



Sp 



c__e 

r 2 



(31) 



The potential given in Eq. ED is repulsive for r < |. For hydrogenic densities 
p 1 ^, the reduced density gradient equals to a ^-independent value of 1.426. 
Near the nuclear cusp, therefore, a pair of electrons on the Is shell provides 
a local repulsive potential which compensates the Coulomb attraction due 
to the nuclear charge. Note that the effective nuclear charge £ for the most 
tightly bound orbital (Is) in a multi-electron atom is smaller than the charge 
of the corresponding nucleus (Z) [40] . As a consequence, the compensation 
is perfect only for one-electron hydrogenic systems. 
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Appendix B 

The effect of approximations to Vt[pA, Pb]{t) on the energies of bound 
states localised in the environment far from subsystem A 

We consider numerical solutions of the Schrodinger equation for one elec- 
tron in the spherically symmetric potential which takes the general from 
given in Eqs. [6J13 The analysed potential corresponds to v^A\fiA',r\ = 0, 
v ext(r) = ~ f> an d v emb !ED [PAiPB] r\ defined in Eq. [7] for p A — > and 
Pb = Pb ■ Such a case represents a local potential around a nucleus with the 
charge Z localised in the environment far from from the investigated subsys- 
tem A and ps comprising entire contribution from doubly occupied Is shell 
centred on this nucleus. Dirac's exchange energy expression (38] is used to 
derive the exchange-correlation component of v^ b , D [pA,PB',r\- The solu- 
tions of one-electron Schrodinger equation for such potential is obtained by 
radial quadrature described in Refs. [531 El] implemented numerically using 
MATLAB2006 environment J55J . Finite difference approximation and the 
matrix representation with 127 radial points is used. The above model sys- 
tem is used to investigate the effect of approximations to v t [pA, Ps]( r ) on the 
lowest energy level. 

Using such approximate potential v t [pA, Ps]( r ), which is finite at nu- 
clear cusp, might lead to appearance of artificially stabilised bound states 
due to improper balance between the nuclear attraction, classical electron- 
electron repulsion, which are both descried exactly and improperly behaving 
Vt[pA, Pb\{t). This imbalance does not cause qualitative problems in the Z=2 
(helium atom) case. No bound states occur if v t [pA, Ps]( r ) is approximated 
by: 



-modelO 



5T nad(TF) [piPB] 



Sp 
-CtfPb 



(32) 

P=PA >0 



derived from the uniform-electron gas expression for v t given in Eq. [TTJ 

For Z = 3 (i.e. Li + ), however, v™ odel ° leads to a bound state of the energy 
-0.2096654 hartree. This value will be used as a reference for analysis of other 
approximations to v t [pA, Pb](^)- 
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The expression derived from the second-order gradient approximation to 
reads: 



sTr d [PA,p B ] 



$PA 



-modell 



5T nad(GEA2)^ 



8p 



(33) 



P=PA >0 

2„ 1 |V7„ 12 



5~ 2/3 1 V PB . 1 |Vp B 

—UtfPb ~ t^t 1 



3 ^ r ^ 72 p B 144 p B 

leads to a further lowering of the energy to -0.611935 hartree. 

In the following part, other potentials will be considered, for which a 
local non-decomposable contribution is added to that given in Eq. [32] near 
the nucleus. 



~model2 i 



r) = -CtfpT (34) 



+ 




i p^. i 



L^pb = L_C1 forr < r T 

^ f- 

for r > r T 



Figure [2] shows that, compared to the decomposable result, the addition of 
a non-decomposable component destabilises the energy but only for such 
small values of r at which sb < 1.6. The maximal destabilisation occurs at 
sb = 1.426 in line with the change of sign of the added term. A similar pic- 
ture emerges from the analysis of the electron density at r = 0. Comparisons 
of Figures [2] and [3] reveals that addition of the non-decomposable term to 
the effective potential affects the orbital energies and electron density in a 
different manner. Even very close to the nucleus this addition reduced elec- 
tron density without affecting the orbital energy noticeably. From the point 
of view of choosing s^ ax determining the range of this additional potential 
it is worthwhile to notice that the limit s = 1.426 should not be exceeded. 
The energy level starts a rapid descend and the on-top density starts to rise 
again at larger values of s%. As far as the lower limit for sb is concerned, 
it should not be smaller than about 0.5 because points at which sb < 0.5 
influence significantly the charge density at the nucleus. 
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Figure 1: Effective potential (Eq. [7]) calculated using local density approx- 
imation for its exchange- and kinetic energy components for: ps = PHe, 

p A — ► 0, and vg t {r) = -2/r. 
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Figure 2: The ground-state energy level in the model system for different 
values of the threshold s# corresponding to the range r T of the | [ pB ' — \ — §^- 

o pB 4 pg 

term in Eq. [3H For s% < 0.376, the corresponding r T does not exists in 
the model system and the results obtained without this additional term are 
shown. 
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Figure 3: The electron density at r = in the model system for different 
values of the threshold s^ corresponding to the range r T of the \ p — \ — §^- 

o pB 4 pg 

term in Eq. [3H For s% < 0.376, the corresponding r T does not exists in 
the model system and the results obtained without this additional term are 
shown. 
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Table 1: Total dipole moment (p in Debye) obtained from the freeze-and- 
thaw calculations using two approximations v t [pA, Pb}^- The target Kohn- 
Sham results are given for comparison. Only dipole moment components 
along the principal axis connecting the subsystems are given. Relative per- 
centage errors of the dipole moments derived from subsystem based calcula- 



tions 



■ A/x 



subsystem _ 



A/J 



An 



Kohn — Sham 



Kohn — Sha >ii 



100%) are given in parentheses. 





freeze- and-thaw 


Kohn-Sham 




Eq.QIl 


Eq.QH 




A 


B 


p 


Li+ 


F" 


5.429 

(9.8) 


5.542 
(7.9) 


6.020 


Li+ 


cr 


5.768 
(15.5) 


5.938 

(13.0) 


6.828 


Li+ 


Br~ 


5.706 
(18.9) 


5.899 

(16.1) 


7.033 


Na+ 


F~ 


7.576 
(1.6) 


7.612 
(1.1) 


7.697 


Na+ 


ci- 


8.509 
(1.2) 


8.563 

(0.5) 


8.608 


Na+ 


Br~ 


8.637 
(1.4) 


8.698 
(0.7) 


8.758 


Be 2+ 


o 2 - 


4.225 
(31.0) 


4.363 

(28.7) 


6.120 


Mg 2+ 


o 2 - 


6.741 
(3.5) 


6.822 
(2.4) 


6.988 


H 2 


H 2 


2.669 

(3.9) 


2.670 
(3.9) 


2.779 


HF 


HF 


3.000 

(2.8) 


3.000 

(2.8) 


3.090 


He 


co 2 


0.014 

(-8.6) 


0.014 

(-8.4) 


0.013 


Ne 


co 2 


0.027 

(-7.6) 


0.027 

(-7.6) 


0.025 
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Table 2: Total dipole moment (/i in Debye) obtained from the freeze- and- 
thaw calculations using two approximations v t \fiA, Pb](t)- The target Kohn- 
Sham results are given for comparison. Only dipole moment components 
along the principal axis connecting the subsystems are given. Relative per- 
centage errors of the dipole moments derived from subsystem based calcula- 
tions 



AfM SU 


bsystem a .K ohn — Sham 
A /.K ohn — Sham 


100%) are 


given ir 


1 parentheses. 








freeze- and-thaw 


Kohn-Sham 








Eq.HU 


Eq.[Tf] 






A 


B 


fi 




Li+ 


H 2 


5.298 


5.353 


5.513 




Li+ 


F 2 


(3.91) 

7.783 


(2.91) 
7.805 


7.807 




Li+ 


co 2 


(0.31) 
9.387 


(0.03) 
9.403 


9.396 




Na+ 


H 2 


(0.10) 
6.512 


(-0.07) 
6.521 


6.527 




H 3 0+ 


Ar 


(0.23) 
2.681 


(0.07) 
2.683 


2.707 




NH+ 


Ar 


(0.96) 
1.661 


(0.89) 
1.665 


1.816 




Be 2 + 


He 


(8.52) 
12.475 


(8.33) 
12.579 


12.986 




Be 2 + 


H 2 


(3.94) 
12.066 


(3.14) 
12.336 


13.569 




Mg 2 + 


He 


(11.08) 
17.302 


(9.09) 
17.313 


17.338 




Mg 2 + 


H 2 


(0.20) 
19.441 
(1.57) 


(0.14) 
19.511 
(1.22) 


19.751 



a for charged systems the dipole moment is ca^ 

origin. 



culated for the cation at the 
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Table 3: The effect of adding non-decomposable contribution to the embed- 
ding potential on the energy on the lowest unoccupied orbital localised in 
the environment (LUEO). For each complex, the used electron density of the 
environment (pb) is the ground-state Kohn-Sham (LDA) electron density of 
the isolated subsystem B. 



subsystem B 


subsystem A 


£lueo 

~TF 
V s — V 


£lueo 
v s = v NDSD 


AeiuEO 


Li+ 


F" 


-2.249 


-2.160 


0.089 


Li+ 


ci- 


-2.669 


-2.568 


0.101 


Li+ 


Br~ 


-2.741 


-2.635 


0.106 


Na+ 


F~ 


-1.908 


-1.874 


0.034 


Na+ 


ci- 


-2.202 


-2.164 


0.038 


Na+ 


Br- 


-2.253 


-2.214 


0.039 


Be 2+ 


o 2 - 


-4.922 


-4.913 


0.009 


Mg 2+ 


o 2 - 


-4.634 


-4.609 


0.025 



Table 4: The effect of adding non-decomposable contribution to the embed- 
ding potential on the energy on the lowest unoccupied orbital localised in 
the environment (LUEO). For each complex, the used electron density of the 
environment (pb) is the ground-state Kohn-Sham (LDA) electron density of 
the isolated subsystem B. 



subsystem B 


subsystem A 


tLUEO 

~TF 
V s — V 


tLUEO 

v s = v NDSD 


^LUEO 


Li+ 


H 2 


-7.533 


-7.332 


0.201 


Li+ 


F 2 


-8.501 


-8.223 


0.279 


Li+ 


co 2 


-9.069 


-8.810 


0.259 


Na+ 


H 2 


-6.295 


-6.240 


0.055 


Be 2 + 


He 


-26.681 


-25.797 


0.884 


Mg 2 + 


He 


-18.376 


-18.078 


0.299 


Be 2 + 

Mg 2 + 


H 2 

H 2 


-21.737 
-16.570 


-21.172 
-16.327 


0.565 
0.242 
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Table 5: The effect of adding non-decomposable contribution to the em- 
bedding potential on the energy on the highest occupied embedded orbital 
(HOEO). For each complex, the used electron density of the environment 
(Pb) is the ground-state Kohn-Sham (LDA) electron density of the isolated 
subsystem B. 

subsystem B subsystem A t HO EO ^hoeo ^hoeo 

v v = v TF v< = v NDSD 



Li+ 


F~ 


-6.716 


-6.635 


0.081 


Li+ 


ci- 


-6.399 


-6.339 


0.060 


Li+ 


Br- 


-6.122 


-6.068 


0.054 


Na+ 


F" 


-5.223 


-5.198 


0.025 


Na+ 


ci- 


-5.397 


-5.379 


0.018 


Na+ 


Br- 


-5.277 


-5.261 


0.016 


Be 2+ 


o 2 - 


-6.724 


-6.690 


0.033 


Mg 2+ 


o 2 - 


-5.135 


-5.108 


0.027 


H 2 O a 


H 2 


-6.799 


-6.799 


0.000 


HF a 


HF 


-9.450 


-9.450 


0.000 


H 2 6 


H 2 


-8.010 


-8.011 


-0.001 


HF 6 


HF 


-10.896 


-10.896 


0.000 



a acceptor of hydrogen bond 
b donor of hydrogen bond 
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Table 6: The effect of adding non-decomposable contribution to the em- 
bedding potential on the energy on the highest occupied embedded orbital 
(HOEO). For each complex, the used electron density of the environment 
(Pb) is the ground-state Kohn-Sham (LDA) electron density of the isolated 
subsystem B. 



subsystem B 


subsystem A 


tHOEO 

~TF 
V s — V 


^HOEO 
v s = V NDSD 


^HOEO 


Li+ 


H 2 


-14.530 


-14.492 


0.037 


Li+ 
Li+ 


F 2 

co 2 


-16.507 
-14.749 


-16.493 
-14.742 


0.014 
0.007 


Na+ 


H 2 


-13.284 


-13.278 


0.006 


Be 2 + 


He 


-37.558 


-37.306 


0.252 


Mg 2 + 


He 


-31.410 


-31.386 


0.024 


Be 2 + 

Mg 2 + 


H 2 

H 2 


-25.208 
-21.178 


-25.001 
-21.131 


0.207 
0.047 
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Table 7: Total dipole moment (/i in Debye) obtained from the freeze- and- 
thaw calculations using different approximations v t [pA, Pb\{t)- The target 
Kohn-Sham results are given in Tables [T] and [2J Only dipole moment com- 
ponents along the principal axis connecting the subsystems are given. Rela- 
tive percentage errors of the dipole moments derived from subsystem based 

' Aii su bsystem A.Kohn — Sham 



calculations 



A/x 



Kohn— Sham 



-100%) are given in parentheses. 



system 


r ,,GEA2 
u t 


„,GEA2 _i_ f n . limit 
u t < J ' u t 


„,GGA97 

u t 


V GGA97 _i_ f . limit 


Li+ 


ci- 


5.365 


5.519 


5.671 


5.829 






(21.4) 


(19.2) 


(16.9) 


(14.6) 


Li+ 


H 2 O a 


_b 


5.188 


5.252 


5.303 






(-) 


(5.9) 


(4.7) 


(3.8) 


Na+ 


ci- 


_b 


8.025 


8.346 


8.400 






(-) 


(6.9) 


(3.0) 


(2.4) 


Na+ 


H 2 O a 


6.378 


6.388 


6.459 


6.470 






(2.3) 


(2.1) 


(1.0) 


(0.9) 


Be 2 + 


o 2 - 


3.962 


4.073 


4.128 


4.245 






(35.3) 


(33.4) 


(32.5) 


(30.6) 


HF 


HF 


2.994 


2.994 


2.995 


2.995 






(3.1) 


(3.1) 


(3.1) 


(3.1) 



a for charged systems the dipole moment is calculated for the cation at the 

origin. b no convergence. 
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